{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Gaussian Processes (GP) with GPy\n",
    "\n",
    "*Deep Bayes summer school, 2018*\n",
    "\n",
    "*A. Zaytsev, Y. Kapushev*\n",
    "\n",
    "\n",
    "In this notebook we are going to use GPy library for GP modeling [SheffieldML github page](https://github.com/SheffieldML/GPy).\n",
    "\n",
    "Why **GPy**?\n",
    "\n",
    "* Specialized library of GP models (regression, classification, GPLVM)\n",
    "* Variety of covariance functions is implemented\n",
    "* There are GP models for large-scale problems\n",
    "* Easy to use"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Run the following line to install GPy library"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "!pip install GPy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "from __future__ import print_function\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import GPy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Current documentation of GPy library can be found [here](http://gpy.readthedocs.org/en/latest/)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Gaussian Process Regression\n",
    "\n",
    "A data set $\\left (X, \\mathbf{y} \\right ) = \\left \\{ (x_i, y_i), x_i \\in \\mathbb{R}^d, y_i \\in \\mathbb{R} \\right \\}_{i = 1}^N$ is given.  \n",
    "\n",
    "Assumption:\n",
    "$$\n",
    "y = f(x) + \\varepsilon,\n",
    "$$\n",
    "where $f(x)$ is a Gaussian Processes and $\\varepsilon \\sim \\mathcal{N}(0, \\sigma_n^2)$ is a Gaussian noise ."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Posterior distribution\n",
    "$$\n",
    "y_* | X, \\mathbf{y}, x_* \\sim \\mathcal{N}(m(x_*), \\sigma(x_*)),\n",
    "$$\n",
    "with predictive mean and variance given by\n",
    "$$\n",
    "m(x_*) = \\mathbf{k}^T \\mathbf{K}_y^{-1} \\mathbf{y} = \\sum_{i = 1}^N \\alpha_i k(x_*, x_i),\n",
    "$$\n",
    "$$\n",
    "\\sigma^2(x_*) = k(x_*, x_*) - \\mathbf{k}^T\\mathbf{K}_y^{-1}\\mathbf{k},\n",
    "$$\n",
    "where\n",
    "$$\n",
    "\\mathbf{k} = \\left ( k(x_*, x_1), \\ldots, k(x_*, x_N) \\right )^T\n",
    "$$\n",
    "$$\n",
    "\\mathbf{K}_y = \\|k(x_i, x_j)\\|_{i, j = 1}^N + \\sigma_n^2 \\mathbf{I}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Building GPR model\n",
    "\n",
    "Lets fit GPR model for function $f(x) = − \\cos(\\pi x) + \\sin(4\\pi x)$ in $[0, 1]$,\n",
    "with noise $y(x) = f(x) + \\epsilon$, $\\epsilon \\sim \\mathcal{N}(0, 0.1)$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "N = 10\n",
    "X = np.linspace(0.05, 0.95, N).reshape(-1, 1)\n",
    "Y = -np.cos(np.pi * X) + np.sin(4 * np.pi * X) + \\\n",
    "    np.random.normal(loc=0.0, scale=0.1, size=(N, 1))\n",
    "plt.figure(figsize=(5, 3))\n",
    "plt.plot(X, Y, '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### 1. Define covariance function\n",
    "\n",
    "The most popular kernel - RBF kernel - has 2 parameters: `variance` and `lengthscale`, $k(x, y) = \\sigma^2 \\exp\\left ( -\\dfrac{\\|x - y\\|^2}{2l^2}\\right )$,\n",
    "where `variance` is $\\sigma^2$, and `lengthscale` - $l$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "input_dim = 1\n",
    "variance = 1\n",
    "lengthscale = 0.2\n",
    "kernel = GPy.kern.RBF(input_dim, variance=variance,\n",
    "                      lengthscale=lengthscale)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### 2. Create GPR model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "model = GPy.models.GPRegression(X, Y, kernel)\n",
    "print(model)\n",
    "model.plot(figsize=(5, 3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Parameters of the covariance function\n",
    "\n",
    "Values of parameters of covariance function can be set like:  `k.lengthscale = 0.1`.\n",
    "\n",
    "Let's change the value of `lengthscale` parameter and see how it changes the covariance function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "k = GPy.kern.RBF(1)\n",
    "theta = np.asarray([0.2, 0.5, 1, 2, 4, 10])\n",
    "figure, axes = plt.subplots(2, 3, figsize=(8, 4))\n",
    "for t, ax in zip(theta, axes.ravel()):\n",
    "    k.lengthscale = t\n",
    "    k.plot(ax=ax)\n",
    "    ax.set_ylim([0, 1])\n",
    "    ax.set_xlim([-4, 4])\n",
    "    ax.legend([t])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Task\n",
    "Try to change parameters to obtain more accurate model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "######## Your code here ########\n",
    "kernel = \n",
    "model = "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "model.Gaussian_noise.variance.fix(0.01)\n",
    "print(model)\n",
    "model.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Tuning parameters of the covariance function\n",
    "\n",
    "The parameters are tuned by maximizing likelihood. To do it just use `optimize()` method of the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "model = GPy.models.GPRegression(X, Y, kernel)\n",
    "model.optimize()\n",
    "print(model)\n",
    "model.plot(figsize=(5, 3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Noise variance\n",
    "\n",
    "Noise variance acts like a regularization in GP models. Larger values of noise variance lead to more smooth model.  \n",
    "Let's check it: try to change noise variance to some large value, to some small value and see the results.\n",
    "\n",
    "Noise variance accessed like this: `model.Gaussian_noise.variance = 1`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "######## Your code here ########\n",
    "model.Gaussian_noise.variance = \n",
    "model.plot(figsize=(5, 3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Now, let's generate more noisy data and try to fit model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "N = 40\n",
    "X = np.linspace(0.05, 0.95, N).reshape(-1, 1)\n",
    "Y = -np.cos(np.pi * X) + np.sin(4 * np.pi * X) + \\\n",
    "    np.random.normal(loc=0.0, scale=0.5, size=(N, 1))\n",
    "\n",
    "kernel = GPy.kern.RBF(1)\n",
    "model = GPy.models.GPRegression(X, Y, kernel)\n",
    "model.optimize()\n",
    "print(model)\n",
    "model.plot(figsize=(5, 3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Now, let's fix noise variance to some small value and fit the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "kernel = GPy.kern.RBF(1)\n",
    "model = GPy.models.GPRegression(X, Y, kernel)\n",
    "model.Gaussian_noise.variance.fix(0.01)\n",
    "model.optimize()\n",
    "model.plot(figsize=(5, 3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Approximate multi-dimensional function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rosenbrock(x):\n",
    "    x = 0.5 * (4 * x - 2)\n",
    "    y = np.sum((1 - x[:, :-1])**2 +\n",
    "                100 * (x[:, 1:] - x[:, :-1]**2)**2, axis=1)\n",
    "    return y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from matplotlib import cm\n",
    "from sklearn.metrics import mean_squared_error\n",
    "\n",
    "def plot_2d_func(func, n_rows=1, n_cols=1, title=None):\n",
    "    grid_size = 100\n",
    "    x_grid = np.meshgrid(np.linspace(0, 1, grid_size), np.linspace(0, 1, grid_size))\n",
    "    x_grid = np.hstack((x_grid[0].reshape(-1, 1), x_grid[1].reshape(-1, 1)))\n",
    "    y = func(x_grid)\n",
    "    fig = plt.figure(figsize=(n_cols * 6, n_rows * 6))\n",
    "    ax = fig.add_subplot(n_rows, n_cols, 1, projection='3d')\n",
    "    ax.plot_surface(x_grid[:, 0].reshape(grid_size, grid_size), x_grid[:, 1].reshape(grid_size, grid_size),\n",
    "                    y.reshape(grid_size, grid_size),\n",
    "                    cmap=cm.jet, rstride=1, cstride=1)\n",
    "    if title is not None:\n",
    "        ax.set_title(title)\n",
    "    return fig"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Here how the function looks like in 2D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plot_2d_func(rosenbrock)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Training set\n",
    "Note that it is 3-dimensional now"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dim = 3\n",
    "x_train = np.random.rand(300, dim)\n",
    "y_train = rosenbrock(x_train).reshape(-1, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task\n",
    "\n",
    "Try to approximate Rosenbrock function using RBF kernel. MSE (mean squared error) should be $<10^{-2}$.\n",
    "**Hint**: if results are not good maybe it is due to bad local minimum. You can do one of the following things:\n",
    "1. Try to use multi-start by calling `model.optimize_restarts(n_restarts)` method of the model.\n",
    "2. Constrain model parameters to some reasonable bounds. You can do it for example as follows:\n",
    "`model.Gaussian_noise.variance.constrain_bounded(0, 1)`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "######## Your code here ########\n",
    "model = "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_test = np.random.rand(3000, dim)\n",
    "y_test = rosenbrock(x_test)\n",
    "y_pr = model.predict(x_test)[0]\n",
    "\n",
    "mse = mean_squared_error(y_test.ravel(), y_pr.ravel())\n",
    "print('\\nMSE: {}'.format(mse))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Covariance functions\n",
    "\n",
    "Short info about covariance function can be printed using `print(k)`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "k = GPy.kern.RBF(1)\n",
    "print(k)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "You can plot the covariance function using `plot()` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "k.plot(figsize=(5, 3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## More \"complex\" functions\n",
    "The most popular covariance function is RBF. However, not all the functions can be modelled using RBF covariance function. For example, approximations of discontinuous functions will suffer from oscillations, approximation of curvy function may suffer from oversmoothing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def heaviside(x):\n",
    "    return np.asfarray(x > 0)\n",
    "\n",
    "\n",
    "def rastrigin(x):\n",
    "    \"\"\"\n",
    "    Parameters\n",
    "    ==========\n",
    "        x : ndarray - 2D array in [0, 1]\n",
    "    \n",
    "    Returns\n",
    "    =======\n",
    "        1D array of values of Rastrigin function\n",
    "    \"\"\"\n",
    "    scale = 8  # 10.24\n",
    "    x = scale * x - scale / 2\n",
    "    y = 10 * x.shape[1] + (x**2).sum(axis=1) - 10 * np.cos(2 * np.pi * x).sum(axis=1)\n",
    "    return y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "fig = plot_2d_func(rastrigin, 1, 2, title='Rastrigin function')\n",
    "\n",
    "x = np.linspace(-1, 1, 100)\n",
    "y = heaviside(x)\n",
    "\n",
    "ax = fig.add_subplot(1, 2, 2)\n",
    "ax.plot(x, y)\n",
    "ax.set_title('Heaviside function')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example of oscillations\n",
    "As you can see there are oscillations in viscinity of discontinuity because we are trying to approximate\n",
    "discontinuous function using infinitily smooth function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "np.random.seed(42)\n",
    "X = np.random.rand(30, 1) * 2 - 1\n",
    "y = heaviside(X)\n",
    "\n",
    "k = GPy.kern.RBF(input_dim=1, variance=1., lengthscale=1.)\n",
    "\n",
    "m = GPy.models.GPRegression(X, y, k)\n",
    "m.optimize()\n",
    "m.plot(figsize=(5, 3))\n",
    "plt.ylim([-0.2, 1.2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example of oversmoothing\n",
    "Actually, the GP model only approximates trend of the function.\n",
    "All the curves are treated as noise.\n",
    "The knowledge about this (in fact there is some repeated structure) should be incorporated into the model via kernel function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(42)\n",
    "\n",
    "X = np.random.rand(300, 2)\n",
    "y = rastrigin(X)\n",
    "\n",
    "k = GPy.kern.RBF(input_dim=2)\n",
    "\n",
    "m = GPy.models.GPRegression(X, y.reshape(-1, 1), k)\n",
    "m.optimize()\n",
    "fig = plot_2d_func(lambda x: m.predict(x)[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Covariance functions in GPy\n",
    "\n",
    "Popular covariance functions: `Exponential`, `Matern32`, `Matern52`, `RatQuad`, `Linear`, `StdPeriodic`. \n",
    "\n",
    "* Exponential:\n",
    "$$\n",
    "k(x, x') = \\sigma^2 \\exp \\left (-\\frac{r}{l} \\right), \\quad r = \\|x - x'\\|\n",
    "$$\n",
    "\n",
    "* Matern32\n",
    "$$\n",
    "k(x, x') = \\sigma^2 \\left (1 + \\sqrt{3}\\frac{r}{l} \\right )\\exp \\left (-\\sqrt{3}\\frac{r}{l} \\right )\n",
    "$$\n",
    "\n",
    "* Matern52\n",
    "$$\n",
    "k(x, x') = \\sigma^2 \\left (1 + \\sqrt{5}\\frac{r}{l} + \\frac{5}{3}\\frac{r^2}{l^2} \\right ) \\exp \\left (-\\sqrt{5}\\frac{r}{l} \\right )\n",
    "$$\n",
    "\n",
    "* RatQuad\n",
    "$$\n",
    "k(x, x') = \\left ( 1 + \\frac{r^2}{2\\alpha l^2}\\right )^{-\\alpha}\n",
    "$$\n",
    "\n",
    "* Linear\n",
    "$$\n",
    "k(x, x') = \\sum_i \\sigma_i^2 x_i x_i'\n",
    "$$\n",
    "\n",
    "* Poly\n",
    "$$\n",
    "k(x, x') = \\sigma^2 (x^T x' + c)^d\n",
    "$$\n",
    "\n",
    "* StdPeriodic\n",
    "$$\n",
    "k(x, x') = \\sigma^2 \\exp\\left ( -2 \\frac{\\sin^2(\\pi r)}{l^2}\\right )\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "covariance_functions = [GPy.kern.Exponential(1), GPy.kern.Matern32(1),\n",
    "                        GPy.kern.RatQuad(1), GPy.kern.Linear(1),\n",
    "                        GPy.kern.Poly(1), GPy.kern.StdPeriodic(1)]\n",
    "figure, axes = plt.subplots(2, 3, figsize=(9, 6))\n",
    "axes = axes.ravel()\n",
    "for i, k in enumerate(covariance_functions):\n",
    "    k.plot(ax=axes[i])\n",
    "    axes[i].set_title(k.name)\n",
    "figure.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Combination of covariance functions\n",
    "\n",
    "* Sum of covariance function is a valid covariance function:\n",
    "\n",
    "$$\n",
    "k(x, x') = k_1(x, x') + k_2(x, x')\n",
    "$$\n",
    "\n",
    "* Product of covariance functions is a valid covariance funciton:\n",
    "$$\n",
    "k(x, x') = k_1(x, x') k_2(x, x')\n",
    "$$\n",
    "\n",
    "### Combinations of covariance functions in GPy\n",
    "\n",
    "In GPy to combine covariance functions you can just use operators `+` and `*`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Let's plot some of the combinations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "covariance_functions = [GPy.kern.Linear(input_dim=1), GPy.kern.StdPeriodic(input_dim=1), GPy.kern.RBF(input_dim=1, lengthscale=1)]\n",
    "operations = {'+': lambda x, y: x + y, '*': lambda x, y: x * y}\n",
    "\n",
    "figure, axes = plt.subplots(len(operations), len(covariance_functions), figsize=(9, 6))\n",
    "\n",
    "import itertools\n",
    "axes = axes.ravel()\n",
    "count = 0\n",
    "for j, base_kernels in enumerate(itertools.combinations(covariance_functions, 2)):\n",
    "    for k, (op_name, op) in enumerate(operations.items()):\n",
    "        kernel = op(base_kernels[0], base_kernels[1])\n",
    "        kernel.plot(ax=axes[count])\n",
    "        axes[count].set_title('{} {} {}'.format(base_kernels[0].name, op_name, base_kernels[1].name),\n",
    "                              fontsize=14)\n",
    "        count += 1\n",
    "figure.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Additive kernels\n",
    "\n",
    "One of the popular approach to model the function of interest is\n",
    "$$\n",
    "f(x) = \\sum_{i=1}^d f_i(x_i) + \\sum_{i < j} f_{ij}(x_i, x_j) + \\ldots\n",
    "$$\n",
    "\n",
    "**Example**: $\\quad f(x_1, x_2) = f_1(x_1) + f_2(x_2)$  \n",
    "To model it using GP use additive kernel $\\quad k(x, y) = k_1(x_1, y_1) + k_2(x_2, y_2)$.\n",
    "\n",
    "More general - add kernels each depending on subset of inputs\n",
    "$$\n",
    "k(x, y) = k_1(x, y) + \\ldots + k_D(x, y),\n",
    "$$\n",
    "where, for example, $k_1(x, x') = k_1(x_1, x_1'), \\; k_2(x, x') = k_2((x_1, x_3), (x_1', x_3'))$, etc.\n",
    "\n",
    "Here is an example of ${\\rm RBF}(x_1) + {\\rm RBF}(x_2)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "k1 = GPy.kern.RBF(1, active_dims=[0])\n",
    "k2 = GPy.kern.RBF(1, active_dims=[1])\n",
    "\n",
    "kernel = k1 + k2\n",
    "\n",
    "x = np.meshgrid(np.linspace(-3, 3, 50), np.linspace(-3, 3, 50))\n",
    "x = np.hstack((x[0].reshape(-1, 1), x[1].reshape(-1, 1)))\n",
    "z = kernel.K(x, np.array([[0, 0]]))\n",
    "\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from matplotlib import cm\n",
    "\n",
    "figure = plt.figure()\n",
    "ax = figure.add_subplot(111, projection='3d')\n",
    "ax.plot_surface(x[:, 0].reshape(50, 50), x[:, 1].reshape(50, 50), z.reshape(50, 50), cmap=cm.jet)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Kernels on arbitrary types of objects\n",
    "\n",
    "Kernels can be defined over all types of data structures: text, images, matrices, graphs, etc. You just need to define similarity between objects.\n",
    "\n",
    "#### Kernels on categorical data\n",
    "\n",
    "* Represent your categorical variable as a by a one-of-k encoding: $\\quad x = (x_1, \\ldots, x_k)$.\n",
    "* Use RBF kernel with `ARD=True`: $\\quad k(x , x') = \\sigma^2 \\prod_{i = 1}^k\\exp{\\left ( -\\dfrac{(x_i - x_i')^2}{\\sigma_i^2} \\right )}$. The lengthscale will now encode whether the rest of the function changes.\n",
    "* Short lengthscales for categorical variables means your model is not sharing any information between data of different categories. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## 2 Sampling from GP\n",
    "\n",
    "So, you have defined some complex kernel.\n",
    "You can plot it to see how it looks and guess what kind of functions it can approximate.\n",
    "Another way to do it is to actually generate random functions using this kernel.\n",
    "\n",
    "GP defines distribution over functions, which is defined by its *mean function* $m(x)$ and *covariance function* $k(x, y)$: for any set $\\mathbf{x}_1, \\ldots, \\mathbf{x}_N \\in \\mathbb{R}^d \\rightarrow$ $\\left (f(\\mathbf{x}_1), \\ldots, f(\\mathbf{x}_N) \\right ) \\sim \\mathcal{N}(\\mathbf{m}, \\mathbf{K})$,\n",
    "where $\\mathcal{m} = (m(\\mathbf{x}_1, \\ldots, \\mathbf{x}_N)$, $\\mathbf{K} = \\|k(\\mathbf{x}_i, \\mathbf{x}_j)\\|_{i,j=1}^N$.\n",
    "\n",
    "Sampling procedure:\n",
    "\n",
    "1. Generate set of points $\\mathbf{x}_1, \\ldots, \\mathbf{x}_N$.\n",
    "2. Calculate mean and covariance matrix $\\mathcal{m} = (m(\\mathbf{x}_1, \\ldots, \\mathbf{x}_N)$, $\\mathbf{K} = \\|k(\\mathbf{x}_i, \\mathbf{x}_j)\\|_{i,j=1}^N$.\n",
    "3. Generate vector from multivariate normal distribution $\\mathcal{N}(\\mathbf{m}, \\mathbf{K})$.\n",
    "\n",
    "Below try to change RBF kernel to some other kernel and see the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "k = GPy.kern.RBF(input_dim=1, lengthscale=0.3)\n",
    "\n",
    "X = np.linspace(0, 5, 500).reshape(-1, 1)\n",
    "\n",
    "mu = np.zeros(500)\n",
    "C = k.K(X, X)\n",
    "\n",
    "Z = np.random.multivariate_normal(mu, C, 3)\n",
    "\n",
    "plt.figure()\n",
    "for i in range(3):\n",
    "    plt.plot(X, Z[i, :])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Task\n",
    "\n",
    "Build a GP model that predicts airline passenger counts on international flights."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "data = np.load('airline.npz')\n",
    "\n",
    "X = data['X']\n",
    "y = data['y']\n",
    "\n",
    "train_indices = list(range(70)) + list(range(90, 129))\n",
    "test_indices = range(70, 90)\n",
    "X_train = X[train_indices]\n",
    "y_train = y[train_indices]\n",
    "\n",
    "X_test = X[test_indices]\n",
    "y_test = y[test_indices]\n",
    "\n",
    "plt.figure(figsize=(5, 3))\n",
    "plt.plot(X_train, y_train, '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You need to obtain something like this\n",
    "\n",
    "<img src=\"airline_result.png\">"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "def plot_model(X, y, model):\n",
    "    x = np.linspace(1948, 1964, 400).reshape(-1, 1)\n",
    "    prediction_mean, prediction_var = model.predict(x)\n",
    "    prediction_std = np.sqrt(prediction_var).ravel()\n",
    "    prediction_mean = prediction_mean.ravel()\n",
    "    \n",
    "    plt.figure(figsize=(5, 3))\n",
    "    plt.plot(X, y, '.', label='Train data')\n",
    "    plt.plot(x, prediction_mean, label='Prediction')\n",
    "    plt.fill_between(x.ravel(), prediction_mean - prediction_std, prediction_mean + prediction_std, alpha=0.3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Let's try RBF kernel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "######## Your code here ########\n",
    "k_rbf = "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see below it doesn't work ;("
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "model = GPy.models.GPRegression(X, y, k_rbf)\n",
    "model.optimize()\n",
    "print(model)\n",
    "plot_model(X_train, y_train, model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "We will try to model this data set using 3 additive components: trend, seasonality and noise.  \n",
    "So, the kernel should be a sum of 3 kernels:  \n",
    "`kernel = kernel_trend + kernel_seasonality + kernel_noise`\n",
    "\n",
    "#### Let's first try to model trend\n",
    "\n",
    "Trend is almost linear with some small nonlinearity, so you can use sum of linear kernel with some other which gives this small nonlinearity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "######## Your code here ########\n",
    "k_trend = "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "model = GPy.models.GPRegression(X, y, k_trend)\n",
    "model.optimize()\n",
    "print(model)\n",
    "plot_model(X_train, y_train, model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Let's model periodicity\n",
    "Just periodic kernel will not work (why?).\n",
    "Try to use product of periodic kernel with some other kernel (or maybe 2 other kernels).\n",
    "Note that the amplitude increases with x."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "######## Your code here ########\n",
    "k_trend = \n",
    "k_seasonal = "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "kernel = k_trend + k_seasonal\n",
    "model = GPy.models.GPRegression(X, y, kernel)\n",
    "model.optimize()\n",
    "print(model)\n",
    "plot_model(X_train, y_train, model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Let's add noise model\n",
    "The dataset is heteroscedastic, i.e. noise variance depends on x: it increases linearly with x.\n",
    "Noise can be modeled using `GPy.kern.White(1)`, but it assumes that noise variance is the same at every x.\n",
    "By what kernel it should be multiplied?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "######## Your code here ########\n",
    "k_trend = \n",
    "k_periodicity = \n",
    "k_noise = "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "kernel = k_trend + k_periodicity + k_noise\n",
    "model = GPy.models.GPRegression(X, y, kernel)\n",
    "model.optimize()\n",
    "print(model)\n",
    "plot_model(X_train, y_train, model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Automatic covariance structure search\n",
    "We can construct kernel is automatic way.\n",
    "Here is our data set (almost the same)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "idx_test = np.where((X[:,0] > 1957))[0]\n",
    "idx_train = np.where((X[:,0] <= 1957))[0]\n",
    "X_train = X[idx_train]\n",
    "y_train = y[idx_train]\n",
    "\n",
    "X_test = X[idx_test]\n",
    "y_test = y[idx_test]\n",
    "\n",
    "plt.figure(figsize=(7, 5))\n",
    "plt.plot(X_train, y_train, '.', color='red');\n",
    "plt.plot(X_test, y_test, '.', color='green');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_model_learned(X, y, train_idx, test_idx, model):\n",
    "    prediction_mean, prediction_var = model.predict(X)\n",
    "    prediction_std = np.sqrt(prediction_var).ravel()\n",
    "    prediction_mean = prediction_mean.ravel()\n",
    "    \n",
    "    plt.figure(figsize=(7, 5))\n",
    "    plt.plot(X, y, '.')\n",
    "    plt.plot(X[train_idx], y[train_idx], '.', color='green')\n",
    "    plt.plot(X, prediction_mean, color='red')\n",
    "    plt.fill_between(X.ravel(), prediction_mean - prediction_std, prediction_mean + prediction_std, alpha=0.3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Expressing Sturcture Through Kernels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example:\n",
    "\n",
    "$$\n",
    "\\underbrace{\\text{RBF}\\times\\text{Lin}}_\\text{increasing trend} + \\underbrace{\\text{RBF}\\times\\text{Per}}_\\text{varying-amplitude periodic} + \\underbrace{\\text{RBF}}_\\text{residual}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Greedy Searching for the Optimum Kernel Combination"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One can wonder: how to automatically search the kernel structure? We can optimize some criteria, which balance between a loss function value and the complexity of the model.\n",
    "Reasinobale candidate for this is BIC-criteria:\n",
    "\n",
    "$$\n",
    "BIC = - 2. \\text{Log-Liklihood} + m \\cdot\\log{n}\n",
    "$$\n",
    "\n",
    "where $n$ sample size and $m$ number of the parameters.\n",
    "\n",
    "However, the procedure of fitting Gaussian Process is quite expensive $O(n^3)$. Hence,  instead of the combinatorial search through all possible combinations, we grow the kernel structure greedy.\n",
    "\n",
    "You can find more details at the https://github.com/jamesrobertlloyd/gp-structure-search. For now, we present toy-example algorithm."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the set of operations:\n",
    "\n",
    "$$\n",
    "\\text{Algebra: } +,\\times\n",
    "$$\n",
    "\n",
    "and the set of basic kernels:\n",
    "\n",
    "$$\n",
    "\\text{Kernels: } \\text{Poly}, \\text{RBF}, \\text{Periodic}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For each level we select extenstion of our current kernel with the lowest BIC. This is an example of the possible kernel grow process (mark notes the lowest BIC at the level):"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src='gp.png'>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task*\n",
    "Implement function that trains a model with given kernel and dataset, calculates and returns BIC\n",
    "The log-lilkelihood of the model can be calculated using `model.log_likelihood()` method,\n",
    "number of parameters of the model you can get via `len(model.param_array)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def train_model_get_bic(X_train, y_train, kernel, num_restarts=1):\n",
    "    '''\n",
    "    Input:\n",
    "        X_train: numpy array of train features, n*d (d>=1)\n",
    "        y_train: numpy array n*1\n",
    "        kernel: GPy object kern\n",
    "        num_restars: number of the restarts of the optimization routine\n",
    "    Output:\n",
    "        bic value\n",
    "    '''\n",
    "    kernel = kernel.copy()\n",
    "    \n",
    "    ######## Your code here ########\n",
    "\n",
    "    return bic "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is a utility function which take list of kernels and operations between them, calculates all product kernels\n",
    "and returns a list of them.\n",
    "After that we need only take sum of the kernels from this list."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def _get_all_product_kernels(op_list, kernel_list):\n",
    "    '''\n",
    "    Find product pairs and calculate them.\n",
    "    For example, if we are given expression:\n",
    "        K = k1 * k2 + k3 * k4 * k5\n",
    "    the function will calculate all the product kernels\n",
    "        k_mul_1 = k1 * k2\n",
    "        k_mul_2 = k3 * k4 * k5\n",
    "    and return list [k_mul_1, k_mul_2].\n",
    "    '''\n",
    "    product_index = np.where(np.array(op_list) == '*')[0]\n",
    "    if len(product_index) == 0:\n",
    "        return kernel_list\n",
    "\n",
    "    product_index = product_index[0]\n",
    "    product_kernel = kernel_list[product_index] * kernel_list[product_index + 1]\n",
    "    \n",
    "    if len(op_list) == product_index + 1:\n",
    "        kernel_list_copy = kernel_list[:product_index] + [product_kernel]\n",
    "        op_list_copy = op_list[:product_index]\n",
    "    else:\n",
    "        kernel_list_copy = kernel_list[:product_index] + [product_kernel] + kernel_list[product_index + 2:]\n",
    "        op_list_copy = op_list[:product_index] + op_list[product_index + 1:]\n",
    "        \n",
    "    return _get_all_product_kernels(op_list_copy, kernel_list_copy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task*\n",
    "\n",
    "This is the main class, you need to implement several methods inside\n",
    "1. method `init_kernel()` - this function constructs initial model, i.e. the model with one kernel. You need just iterate through the list of base kernels and choose the best one according to BIC\n",
    "2. method `grow_level()` - this function adds new level. You need to iterate through all base kernels and all operations,\n",
    "apply each operation to the previously constructed kernel and each base kernel (use method `_make_kernel()` for this) and then choose the best one according to BIC."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class GreedyKernel:\n",
    "    '''\n",
    "    Class for greedy growing kernel structure\n",
    "    '''\n",
    "    def __init__(self, algebra, base_kernels):\n",
    "        self.algebra = algebra\n",
    "        self.base_kernels = base_kernels\n",
    "        self.kernel = None\n",
    "        self.kernel_list = []\n",
    "        self.op_list = []\n",
    "        self.str_kernel = None\n",
    "    \n",
    "    def _make_kernel(self, op_list, kernel_list):\n",
    "        '''\n",
    "        Sumation in kernel experssion\n",
    "        '''\n",
    "        kernels_to_sum = _get_all_product_kernels(op_list, kernel_list)\n",
    "        new_kernel = kernels_to_sum[0]\n",
    "        for k in kernels_to_sum[1:]:\n",
    "            new_kernel = new_kernel + k\n",
    "        return new_kernel\n",
    "    \n",
    "    def init_kernel(self, X_train, y_train):\n",
    "        '''\n",
    "        Initialization of first kernel\n",
    "        '''\n",
    "        best_kernel = None\n",
    "        \n",
    "        ###### Your code here ######\n",
    "\n",
    "        # You need just iterate through the list of base kernels and choose the best one according to BIC\n",
    "        # save the kernel in `best_kernel` variable\n",
    "        \n",
    "        # base kernels are given by self.base_kernels --- list of kernel objects\n",
    "        \n",
    "        \n",
    "        \n",
    "        \n",
    "        ############################\n",
    "        \n",
    "        assert best_kernel is not None\n",
    "        \n",
    "        self.kernel_list.append(best_kernel)\n",
    "        self.str_kernel = str(best_kernel.name)\n",
    "        \n",
    "    def grow_level(self, X_train, y_train):\n",
    "        '''\n",
    "        Select optimal extension of current kernel \n",
    "        '''\n",
    "        \n",
    "        best_kernel = None  # should be kernel object\n",
    "        best_op = None  # should be operation name, i.e. \"+\" or \"*\"\n",
    "        \n",
    "        ###### Your code here ######\n",
    "        \n",
    "        # You need to iterate through all base kernels and all operations,\n",
    "        # apply each operation to the previously constructed kernel and each base kernel\n",
    "        # (use method `_make_kernel()` for this) and then choose the best one according to BIC.\n",
    "        \n",
    "        # base kernels are given by self.base_kernels --- list of kernel objects\n",
    "        # operations are given by self.algebra --- dictionary:\n",
    "        #                                              {\"+\": lambda x, y: x + y\n",
    "        #                                               \"*\": lambda x, y: x * y}\n",
    "\n",
    "        # best_kernel - kernel object, store in this variable the best found kernel\n",
    "        # best_op - '+' or '*', store in this variable the best found operation\n",
    "        \n",
    "        \n",
    "        \n",
    "        \n",
    "        ############################\n",
    "\n",
    "        assert best_kernel is not None\n",
    "        assert best_op is not None\n",
    "        \n",
    "        self.kernel_list.append(best_kernel)\n",
    "        self.op_list.append(best_op)\n",
    "        \n",
    "        new_kernel = self._make_kernel(self.op_list, self.kernel_list)\n",
    "        str_new_kernel = '{} {} {}'.format(self.str_kernel, best_op, best_kernel.name)\n",
    "        \n",
    "        return new_kernel, str_new_kernel\n",
    "    \n",
    "    def grow_tree(self, X_train, y_train, max_depth):\n",
    "        '''\n",
    "        Greedy kernel growing\n",
    "        '''\n",
    "        if self.kernel == None:\n",
    "            self.init_kernel(X_train, y_train)\n",
    "            \n",
    "        for i in range(max_depth):\n",
    "            self.kernel, self.str_kernel = self.grow_level(X_train, y_train)\n",
    "            print(self.str_kernel)\n",
    "            \n",
    "    def fit_model(self, X_train, y_train, kernel, num_restarts=1):\n",
    "        model = GPy.models.GPRegression(X_train, y_train, kernel)\n",
    "        model.optimize_restarts(num_restarts, verbose=False)\n",
    "        return model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let us define the algebra and list of base kernels.\n",
    "To make learning process more robust we constrain some parameters of the kernels to lie within\n",
    "some reasonable intervals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# operations under kernels:\n",
    "\n",
    "algebra = {'+': lambda x, y: x + y,\n",
    "           '*': lambda x, y: x * y\n",
    "          }\n",
    "\n",
    "# basic kernels list:\n",
    "poly_kern = GPy.kern.Poly(input_dim=1, order=1)\n",
    "\n",
    "periodic_kern = GPy.kern.StdPeriodic(input_dim=1)\n",
    "periodic_kern.period.constrain_bounded(1e-2, 1e1)\n",
    "periodic_kern.lengthscale.constrain_bounded(1e-2, 1e1)\n",
    "\n",
    "rbf_kern = GPy.kern.RBF(input_dim=1)\n",
    "rbf_kern.lengthscale.constrain_bounded(1e-2, 1e1)\n",
    "\n",
    "bias_kern = GPy.kern.Bias(1)\n",
    "\n",
    "kernels_list = [poly_kern, periodic_kern, rbf_kern]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's train the model.\n",
    "You should obtain something which is more accurate than the trend model ;)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "GK = GreedyKernel(algebra, kernels_list)\n",
    "GK.grow_tree(X_train, y_train, 4)\n",
    "model = GK.fit_model(X_train, y_train, GK.kernel)\n",
    "plot_model_learned(X, y, idx_train, idx_test, model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bonus Task\n",
    "Try to approximate rastrigin function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plot_2d_func(rastrigin)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Training set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.rand(42)\n",
    "x_train = np.random.rand(200, 2)\n",
    "y_train = rastrigin(x_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Hint: you can constrain parameters of the covariance functions, for example\n",
    "`model.std_periodic.period.constrain_bounded(0, 0.2)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "######## Your code here ########\n",
    "model = "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(model)\n",
    "\n",
    "x_test = np.random.rand(1000, 2)\n",
    "y_test = rastrigin(x_test)\n",
    "y_pr = model.predict(x_test)[0]\n",
    "\n",
    "mse = mean_squared_error(y_test.ravel(), y_pr.ravel())\n",
    "print('MSE: {}'.format(mse))\n",
    "\n",
    "fig = plot_2d_func(lambda x: model.predict(x)[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Appendix: Gaussian Process Classification"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Classification\n",
    "\n",
    "A data set $\\left (X, \\mathbf{y} \\right ) = \\left \\{ (x_i, y_i), x_i \\in \\mathbb{R}^d, y_i \\in \\{+1, -1\\} \\right \\}_{i = 1}^N$ is given.  \n",
    "\n",
    "Assumption:\n",
    "$$\n",
    "p(y = +1 \\; | \\; x) = \\sigma(f(x)) = \\pi(x),\n",
    "$$\n",
    "where latent function $f(x)$ is a Gaussian Processes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "We need to produce a probabilistic prediction\n",
    "$$\n",
    "\\pi_* = p(y_* \\; | \\; X, \\mathbf{y}, x_*) = \\int \\sigma(f_*) p(f_* \\; | \\; X, \\mathbf{y}, x_*) df_*,\n",
    "$$\n",
    "$$\n",
    "p(f_* \\; | \\; X, \\mathbf{y}, x_*) = \\int p(f_* \\; | \\; X, x_*, \\mathbf{f}) p(\\mathbf{f} \\; | \\; X, \\mathbf{y}) d\\mathbf{f},\n",
    "$$\n",
    "where $p(\\mathbf{f} \\; |\\; X, \\mathbf{y}) = \\dfrac{p(\\mathbf{y} | X, \\mathbf{f}) p(\\mathbf{f} | X)}{p(\\mathbf{y} | X)}$ is the posterior over the latent variables.\n",
    "\n",
    "Both integrals are intractable.\n",
    "\n",
    "Use approximation technique like Laplace approximation or Expectation Propagation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "from matplotlib import cm\n",
    "\n",
    "def cylinder(x):\n",
    "    y = (1 / 7.0 - (x[:, 0] - 0.5)**2 - (x[:, 1] - 0.5)**2) > 0\n",
    "    return y\n",
    "\n",
    "np.random.seed(42)\n",
    "X = np.random.rand(40, 2)\n",
    "y = cylinder(X)\n",
    "\n",
    "x_grid = np.meshgrid(np.linspace(0, 1, 100), np.linspace(0, 1, 100))\n",
    "y_grid = cylinder(np.hstack((x_grid[0].reshape(-1, 1), x_grid[1].reshape(-1, 1)))).reshape(x_grid[0].shape)\n",
    "\n",
    "positive_idx = y == 1\n",
    "plt.figure(figsize=(5, 3))\n",
    "plt.plot(X[positive_idx, 0], X[positive_idx, 1], '.', markersize=10, label='Positive')\n",
    "plt.plot(X[~positive_idx, 0], X[~positive_idx, 1], '.', markersize=10, label='Negative')\n",
    "im = plt.contour(x_grid[0], x_grid[1], y_grid, 10, cmap=cm.hot)\n",
    "plt.colorbar(im)\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "kernel = GPy.kern.RBF(2, variance=1., lengthscale=0.2, ARD=True)\n",
    "\n",
    "model = GPy.models.GPClassification(X, y.reshape(-1, 1), kernel=kernel)\n",
    "model.optimize()\n",
    "print(model)\n",
    "\n",
    "\n",
    "def plot_model_2d(model):\n",
    "\n",
    "    model.plot(levels=40, resolution=80, plot_data=False, figsize=(5, 3))\n",
    "    plt.plot(X[positive_idx, 0], X[positive_idx, 1], '.', markersize=10, label='Positive')\n",
    "    plt.plot(X[~positive_idx, 0], X[~positive_idx, 1], '.', markersize=10, label='Negative')\n",
    "    plt.legend()\n",
    "    plt.show()\n",
    "    \n",
    "plot_model_2d(model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Let's change lengthscale to some small value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "model.rbf.lengthscale = [0.05, 0.05]\n",
    "plot_model_2d(model)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
